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At a quantum critical point, bipartite entanglement entropies have universal quantities which are 
subleading to the ubiquitous area law. For Renyi entropies, these terms are known to be similar to 
the von Neumann entropy, while being much more amenable to numerical and even experimental 
measurement. We show here that when calculating universal properties of Renyi entropies, it is 
important to account for unusual corrections to scaling that arise from relevant local operators 
present at the conical singularity in the multi-sheeted Riemann surface. These corrections grow in 
importance with increasing Renyi index. We present studies of Renyi correlation functions in the 
l-fl transverse-field Ising model (TFIM) using conformal field theory, mapping to free fermions, and 
series expansions, and the logarithmic entropy singularity at a corner in 2 + 1 for both free bosonic 
field theory and the TFIM, using numerical linked cluster expansions. In all numerical studies, 
accurate results are only obtained when unusual corrections to scaling are taken into account. 

In the worst case, an analysis ignoring these corrections can get qualitatively incorrect answers, 
such as predicting a decrease in critical exponents with the Renyi index, when they are actually 
increasing. We discuss a two-step extrapolation procedure that can be used to account for the 
unusual corrections to scaling. 


I. INTRODUCTION 

Quantum entanglement in the ground state of a many- 
body system contains universal signatures of a vari¬ 
ety of low-energy and long-lengthscale physical phenom¬ 
ena. These include Goldstone modes due to spontaneous 
breaking of a continuous symmetry [IHH; topological or¬ 
der m-, Fermi surfaces mm-, quantum criticality nni- 
I13j : and more. If a system is divided into two spatial 
regions A and B, the entanglement in a pure state can 
be characterized by the von Neumann entropy associated 
with the reduced density matrix of either subsystem, 

Si =-TrpA\^PA- (1) 

However, numerically computing the von Neumann en¬ 
tropy is a significant challenge, requiring an explicit rep¬ 
resentation of the reduced density matrix or the ground 
state wavefunction. In continuum field theories, the 
von Neumann entropy can be computed from the Renyi 
entropies, 

5a = :p^ln(TrpS), (2) 

I — a 

by introducing a identical replicas of the system, then 
formally taking the limit a —>■ 1 m- Renyi entropies for 
integer values of a > 2 can be calculated by a variety of 
methods employing this replica-trick, however, the limit 
a —>■ I is difficult to obtain from numerical data for finite, 
discrete a, which can often be most accurate for small a. 


For this reason, integer a > 1 Renyi entropies are of¬ 
ten studied in quantum many-body systems on their own 
merit. A variety of recent numerical techniques have been 
developed for this task. One way to calculate them in 
Quantum Monte Carlo (QMC) simulations is by evaluat¬ 
ing the expectation value of a swap operator between the 
different replicas m- Alternatively, these integer Renyi 
entropies can be expressed as a ratio of two partition 
functions [16] and evaluated by stochastically sampling 
an extended ensemble where one switches between the 
configurations defining the two partition functions m 
Such sampling techniques make the Renyi entropies avail¬ 
able in other numerical approximations of the wave func¬ 
tion, such as tensor network states [TS] . They can also be 
calculated by series expansions in various coupling con¬ 
stants of the model diiiMa- Most importantly, there 
are a number of proposals [23H25] to measure the Renyi 
entropy experimentally. 

This manuscript considers the universal properties as¬ 
sociated with the Renyi entropies and associated Renyi 
correlation functions in a critical system. While Renyi 
entropies do not satisfy strong subadditivity a priori, and 
thus may not provide as rigorous an information mea¬ 
sure of entanglement as Eq. ([^, explicit calculations have 
found that their singular behavior is very analogous to 
the von Neumann entropy in typical many-body systems. 
Universal terms associated with Goldstone modes, topo¬ 
logical order, Fermi surfaces and quantum critical points 
often depend in a simple parametric way on the Renyi 
index a. Hence, computing universal terms in any one 
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Renyi entropy can be sufficient to determine the univer¬ 
sality class of the system. If there is an entanglement 
property that exhibits monotonicity under renormaliza¬ 
tion [26j . one would expect it to remain valid for the 
Renyi entropies as well. 

However, there is a subtlety in obtaining universal, 
asymptotic Renyi properties, which is the main focus of 
this work. Renyi entropies, and associated correlation 
functions, suffer from “unusual” corrections to scaling 
that are absent for von Neumann entropies [57]. These 
corrections become more and more severe as the Renyi in¬ 
dex increases and cause a slow convergence to the asymp¬ 
totic limit. Rapid convergence in the von Neumann en¬ 
tropy and slower convergence for higher Renyi cases were 
previously noted for topological entanglement entropy, 
where they were found to be related to non-universal per¬ 
turbative terms [55]. Here, we discuss a two step extrap¬ 
olation process needed to obtain the asymptotic universal 
critical Renyi properties. 

In light of our work, there may be a need to revisit 
previous calculations of Renyi entropies and their uni¬ 
versal terms in past computational studies. For exam¬ 
ple, in this paper we re-examine recent numerical studies 
of 0{N) models in 2 -|- 1 dimension, where it was found 
that the coefficients of the logarithmic singularities asso¬ 
ciated with a corner approximately scales with N [5M55] . 
These studies also found signihcant differences between 
the entropy of N free fields and a single copy of the inter¬ 
acting 0{N) theory (differences of order of greater than 
10 percent). In light of a recent conjecture that the cor¬ 
ner entanglement mimics a central charge Ct na , one 
might expect a much smaller difference. Here, we find 
that some of the observed difference could be due to the 
unusual corrections to scaling. Our work also suggests 
that one should be especially careful in computing uni¬ 
versal properties at large values of the Renyi index a. 

Following a discussion of how unusual corrections to 
Renyi properties arise in Section |Tl| we illustrate their 
importance in Section in by studying Renyi correlation 
functions in (space-l-time) dimension D = 1-1-1, for which 
there is a well-developed theory. In Section jlV] we apply 
these insights in H = 2 -|- 1 to universal entanglement en¬ 
tropy contributions from sharp corners. Accounting for 
the unusually strong Renyi finite-size corrections signif¬ 
icantly improves our numerical results for the universal 
corner term for the free boson theory; we apply the same 
procedure to improve our estimates of the corner term 
for the critical Ising model. 


II. UNUSUAL CORRECTIONS TO SCALING 


Corrections to critical scaling are well known for bulk 
quantities and are due to the presence of irrelevant oper¬ 
ators [S.lj . When approaching a critical point, a thermo¬ 
dynamic quantity such as the order parameter suscepti¬ 
bility X does not behave as a single power-law, but as a 
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FIG. 1. A schematic view of the conical singularity for the 
a = 2 Renyi partition function TZ 2 - Dashed lines represent the 
periodic boundaries in imaginary time. All lattice sites retain 
a coordination number of four, however the re-arrangement 
of periodic boundaries in imaginary time near the conical sin¬ 
gularity (red) push that local region off-critical. 


sum of power laws. For small non-zero t (x T — Tc, 

xit)=At-^il + Bt^' +...), (3) 


or for a finite system of size L at the critical point, 

xiL)=AL^'il + C/L^ + ...). (4) 


The subleading terms above lead to violations of simple 
scaling and data collapse, such as the failure of Binder 
ratios to cross precisely at a critical point. The sublead¬ 
ing terms can be seen to directly affect the calculation of 
exponents, by taking a logarithm of Eq. (j^ and differen¬ 
tiating with respect to t for small t, 


d\ogX 

dlogt 




(5) 


Thus the effective exponent depends on t and, by scaling 
theory, on L, 

7eff = 7 + S"/ = 7 + C/L^ (6) 

Note that the same is true for the coefficient of a loga¬ 
rithmic singularity. 

When studying lattice theories in the bulk, such as the 
3-dimensional classical Ising model, including or ignor¬ 
ing such conventional scaling corrections in the analysis 
of finite-size data typically makes a difference of less than 
one percent in estimates of critical exponents such as 7 . 
For this reason, analyses of Monte Carlo simulation or se¬ 
ries expansion results often ignore subleading corrections 
to finite-size scaling. However, as first discussed by Cardy 
and Calabrese [27] in the context of 1 -I- 1-dimensional 
conformal field theories (CFTs), unusual corrections to 
scaling arise in the calculation of Renyi entropies. This 
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can be understood in the path-integral picture, where the 
entropy is proportional to the partition function of a a- 
sheeted Riemann surface 'R-a. [131133] (Fig. [^, which at 
finite-temperature involves a branch cut located at nf3 for 
n = 1,2, ■ ■ ■ a between each of the a replicas defined on 
region B. A conical singularity exists on the boundary 
between regions A and B where this branch cut intersects 
the a/3-periodic imaginary time structure of the path in¬ 
tegral in region A. Local operators, which might even 
be relevant in the bulk, are introduced at this conical 
singularity in the multi-sheeted Riemann surface. Even 
though their contribution is only local, these operators 
yield subleading corrections to scaling with a power that 
nevertheless diminishes with the Renyi index a and ulti¬ 
mately goes to zero as a goes to infinity, strongly altering 
the scaling behavior. 

The unusual correction to scaling exponents in 
D = 1 \ are well known m and have been clearly ob¬ 

served in numerical and analytical lattice calculations 
of spin chains [55H55] . For the case of the D = 1 -|- 1 
transverse-held Ising model, which we use as an example 
in the next section, the most relevant operator that does 
not break the Z 2 symmetry is the energy operator with 
a scaling dimension Xe = ^- For the Renyi correlation 
functions (defined in the next section) it leads to a finite- 
size correction to scaling like that of Eq. Q but with an 
exponent. 


9 = xeIol, 


(7) 


where are the Pauli operators. We focus on the 

critical point h/ J = 1, or its vicinity. The effect of un¬ 
usual corrections to scaling has been well studied in this 
model for the Renyi entropies |3H] and related quantities 
(321 US]. Here, we instead examine the scaling of esti¬ 
mators weighted by the Renyi moments. One advantage 
of looking at these Renyi correlations is that one can 
observe more directly how the conical singularity affects 
correlations, not just the partition function (the entropy). 
As we will demonstrate, these correlation can easily be 
computed with conformal held theory techniques. 

To dehne the Renyi correlators, we divide the system 
into two halves A and B and calculate the expectation 
value of operators (or pairs of operators) in the Renyi 
multi-sheeted geometry TZa- The expectation value of an 
observable O in subsystem A is dehned as. 


(O) 


Tr (p^)“ O 
Tr (/3 a )“ 


( 10 ) 


where pA is the reduced density matrix for subsystem A 
and a is the Renyi index. Note that setting a = 1 gives 
us normal bulk expectation values, which can also be 
obtained without partitioning of the system, providing a 
non-trivial check on our computational procedures. 


A. Renyi correlators in conformal field theory 


which goes to zero for large a. 

While there is no theory for the unusual corrections 
to scaling in Zl = 2 -f 1, it is reasonable to expect that 
they are again dominated by the presence of the energy 
operator at the conical singularity. The scaling dimension 
of the energy operator is related to the correlation 
length exponent e by the relation m, 

xe = d + 1 — l/v. (8) 


We consider the simplest D = 1 \ setup, which is a 

subsystem of length ^ in an infinite system. As usual for 
such computations we use the replica trick: a is hrst as¬ 
sumed to be an integer, but the formulae are expected to 
hold for any a > 0. As explained in e.g. Refs. [MllM], the 
entropy is - up to some normalization - nothing but the 
partition function of a a-sheeted Riemann surface TZa- 
The Riemann surface may be mapped onto the complex 
plane by the conformal transformation, 


In particular, = 1 for the D = 2 + 1 free boson theory, 
and approximately 1.41 for the 2-|-l Ising theory |3I]. In 
D = 2 + 1 there is no particular reason why the correction 
should be given by the simple formula XeIol, as in 13 = 
1-1-1. In Section [rv] we will see that while the exponent 
decreases with a, the numerical results suggest that the 
dependence on a is slightly different. 


III. RENYI CORRELATORS IN D = 1 -t 1 

To illustrate a case where failing to account for the 
unusual corrections to scaling can lead to qualitative 
changes in finite-size scaling, consider the transverse-field 
Ising model (TFIM) in spatial dimension d = 1. The 
Hamiltonian is, 

H = J (9) 

i i 


Here, w is a complex coordinate on the Riemann surface 
TZa, while z is a complex coordinate in the complex plane 
C. Subsystem A corresponds to w real in the interval 
[ 0 ,£]. 

To obtain a Renyi correlator in the ground state, what 
we need is to compute the correlation on TZa- Let us do 
that for some two point function of a primary operator. 
We use the two point function of a primary operator with 
scaling dimension x, together with the transformation 
law of a primary operator |44j . We obtain, 

(4>(u;)4>(u;'))^^ = 

Here z (resp. z') has to be understood as z{w) (resp. 


dz dz' 


dw dw' (z — z')^ 


( 12 ) 
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z{w')). Using (111 this becomes, 




£ 

( WW' \2 + 2 / 

y (£—w)(£—w') J 

avow' 

( w \°‘ 

/ u;' \ “ 


1 ^ 

1 

\l-w' ) 


(13) 

Eq. ^ is the main result of this section. The Renyi 
correlations are translational invariant for a = 1, as they 
should be. This is also the case far from the boundaries 
w = 0,i for any a. Similar to the entanglement entropy 
[H, various generalizations to finite systems and finite 
temperatures are straightforward [IK] . 

In the following we are going to check our results in the 
Ising chain. The most natural example is the case when 
one of the points lies at the end of the interval wi = S {S 
is a UV cutoff of the order of a lattice spacing; it ensures 
that the correlations stay finite), and the other at some 
distance r from it: 


($((5)$(r))K„ 


£ 

( \ 

l/2-|-l/2a 

J 

aSr , 

\ 

5,1/a 

tJ 


1/a 


(14) 


(1 


\ X j Ol — X 

~ j) 

j.xjcx-\-x 


(^Aa + 0{r “)),(15) 


where is some a-dependent UV cutoff. The correction 
terms 0(r“^/“) lead to the correction to scaling exponent 
1/a in Eq. ([^. We are going to test formula (15) in 
the following subsections. Note also, we recover that the 
scaling dimension of an operator near the cone is modified 
to X —>■ X/a (Eq. 0)> as pointed out in Ref. m 


B. A numerical check using free fermions 


It is well known that the Hamiltonian 0 can be 
mapped onto a system of free fermions through a Jordan- 
Wigner transformation. 


- 1 , 

2 cr+ = CTj + iaj = 2cJ exp 



(16) 

(17) 


To perform the identihcation one of the easiest way is to 
use the correlation matrix, 


K = 


C 


D 


D* l-C* 


(18) 


where M* (resp. M*) denotes the transpose (resp. con¬ 
jugate transpose) of any matrix M. The f x € matrices 
C and D encode the two types of correlators: 


C — ((cJcj))i<ij<.f, 
D = ((cJcJ-))i<,,j<^ 


(19) 

( 20 ) 


The eigenvalues of K are simply related to the single 
particle eigenenergies of the entanglement Hamiltonian 
|4H] , and the knowledge of the former allows to determine 
the latter. Said differently, given a subsystem correlation 
matrix it is easy to determine the corresponding entan¬ 
glement Hamiltonian, and then to compute the entropy. 

Here we want to compute the Renyi correlations, as 
defined in Eq. (10). Since the RDM is the exponential 


of a quadratic form, so is any power of it. The new en¬ 
tanglement Hamiltonian corresponding to the new RDM 
pJ^/TrpJ^ is then simply aHj^. One can then find the 
modified Renyi correlation matrix Ka that would have 
uHa as the entanglement Hamiltonian. The details are 
explained in Appendix In compact form the result 
may be written. 


= 


l + (A-i-l)^ 


( 21 ) 


Note that in absence of pairing terms c\ the same for¬ 
mula holds with K and replaced by C and Ca- Us¬ 
ing (21), one can access all Renyi correlations between 

the spin correlations in the z 


fermions. 


(4o)c 


(44 


basis being recovered by using the Jordan Wigner trans¬ 
formation (16 17), and appropriate treatment of the re¬ 
sulting Jordan-Wigner string (see e.g. Ref. fiKIl . For the 
spin-spin correlation at distance r we obtain. 


cr,- cr. 


+r 


)„ — {{ai+rn-lbi+n) c 


( 22 ) 


with aj = c] — Cj and bj = Cj + cj. 


We now use (22) to numerically compute the correla¬ 
tion, 


and then diagonalized by a Bogoliubov transformation. 
The correlations in the ground state may then be ob¬ 
tained analytically [iH] . 

It is also well known how to compute the entangle¬ 
ment entropy for such quadratic fermion systems mm- 
The method relies on the following observation. Since 
Wick’s theorem holds for any fermionic correlation in 
subsystem A, the reduced density matrix (RDM) itself 
is the exponential of a quadratic fermion Hamiltonian 
Ha, Pa = /Tv . Finding the precise form of 
Ha can then be done by demanding that it reproduce all 
correlation functions in subsystem A, (O) = Tr(p^O). 


a(r) = 


(23) 


in the ground state of the chain 0 at the critical point 
h/J=l. We consider several values of a and several dis¬ 
tances 0 < r < £ in a subsystem of length £ = 4096. As 
can be seen in Fig. the agreement with the prediction 
(15) with Ising/Onsager exponent a: = 1/8 is very good. 
Note however that finite-size effects increase significantly 
with a; indeed the formula is expected to be valid only in 
the regime ^ 1, which for moderately large values 
of a already requires immense system sizes. Another in¬ 
teresting feature is that the correlation near r = £ blows 
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FIG. 2. Renyi spin-spin correlation Ca{r) = (croa^)^ in the 
Ising chain for several values of a = 0.5,1,2, 3. The full sys¬ 
tem is infinite, and the subsystem size is i = 4096. The dots 
are the numerical data, while the lines are the CFT result. 
The UV cutoff Ac in (151 is adjusted such that the analytical 
result matches the numerical one at r = ^/2. 


up again at r = £ for a > 1. One possible interpretation 
would be that since the scaling dimension x near the cone 
is modified to x/a, the spin configurations near the (fic¬ 
titious) endpoints at r = 0 ,.^ are more constrained, so 
more likely to match. Since translational invariance is 
broken, the effect appears similar to that of a boundary 
with partial reflection. 

Of course, the fact that we are able to test our re¬ 
sults with good precision relies heavily on the free fermion 
structure. For more general system the computation of 
Renyi correlations is more difficult, and the accessible 
system sizes are much smaller. To illustrate this we look 
at similar correlations with a different method below. 


C. Correlation Function and Structure Factors 


Let us calculate correlation functions Ca(r) directly in 
the high-field (disordered) phase of Eq. ([^ by a series ex¬ 
pansion in J/h [3^. We define Ca{r) via Eq. (23), where 
0 is a boundary site of A directly at the interface with B 
and r refers to any other site in subsystem A. Eor corre¬ 
lation function at r in Eq. (23), the series coefficients will 
be zero below order r. One needs at least r powers of the 
perturbation to get a non-zero correlation between site 0 
and site r. This means as one goes further away in dis¬ 
tance, the series become effectively shorter and shorter. 
Eor this reason, analyzing real-space correlations is not 
very useful. Series extrapolation methods, which are nec¬ 
essary for studying the critical region, simply cannot be 
applied. Instead, we need to focus on structure factors 


or correlation sums defined as. 


5„=^C„(r). 


(24) 


It is useful but not necessary to multiply all r > 0 terms 
in the sum in Eq. (24) by a factor of 2, as that would 


be the equal-time structure factor of the ID system in 
the absence of a partition, where one has two neighbors 
at each such distance. This is what we have done in the 
analysis discussed below. 

We describe the series analysis using the method of 
differential approximants [501 EH in Appendix [b] There, 
we define 0 = M-\-N-\-J+l as the order of the series 
used in the estimation of the critical exponent (M, iV, 
and J are individually the orders of three polynomials 
defining the extrapolation scheme). In other words, O is 
the maximum power used in determining the polynomial 
coefficients. The order of the series expansion is related 
to a length scale [50], thus effective exponents can vary 
with the order as. 


7eff = 7 + C/0®. 


(25) 


Again, theory predicts that the scaling dimension of the 
bulk spin is given by, x^uik = 1/8. The scaling dimension 
of the boundary spin with Renyi index a is given by. 


^Ol — 


^bulk 


a 


1 

8 a 


(26) 


Thus the correlation function should decay for large r as, 

C^{r) = (27) 

Summing over r, the structure factor Sa should diverge 
on approach to the critical point with an exponent. 


Tct (1 ^bulk XQ^')y. 


(28) 


Since, zz = 1 for the model, the exponent simplifies to, 

7a = (1 - Xbulk - aza) = (1 - ^ - ^)- (29) 

8 8 a 

Figure a) shows the estimated structure factor expo¬ 
nent obtained from order (J//i)^^ series expansions with 
various choices for the order N of one of the polyno¬ 
mials used in the fit. Note that the exponent 7 is con¬ 
verged very well to the theoretical value for a = I (which 
gives the bulk exponent), but the convergence is poor for 
larger a. In fact, while theory predicts that the expo¬ 
nent 7 should increase with a, the series extrapolation 
estimate appears to systematically decrease. 

Figure [^b) shows plots of the estimated exponent val¬ 
ues for different orders of the series. The x axis has been 
scaled differently for different a values. The rescaling 
pushes the data for larger a values more and more to the 
right. Thus for larger a one needs to extrapolate farther 
to reach the asymptotic value. The solid squares are the 
asymptotic theoretical values and the dashed lines are a 
guide to the eye showing the results are clearly consistent 
with theory, only when the unusual corrections to scaling 
are taken into account. 
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FIG. 3. (a) Structure factor exponents estimated with the full length series. The A'^-axis labels different approximates. The 
solid circles are the estimated values (with dashed lines as a guide to the eye). The solid lines are the theoretically predicted 
values, (b) Estimated exponents for different orders O of the series expansion. The squares are the asymptotic theoretical 
values. Dashed lines are guides to the eye, showing that a rough extrapolation of data is consistent with theory. 


D. Magnetization in the ordered phase 

We now calculate series expansions for the Renyi mag¬ 
netization at the boundary site of A by a series expansion 
in h/J (a low-field expansion in the ordered phase). This 
boundary magnetization is defined as, 

AfO = (ag)„. (30) 

Once again, setting a = 1 just gives us the usual bulk 
magnetization that is well known to vanish at the critical 
point with exponent /3 = 1/8. 

When extrapolating, we do not allow any P^{x) terms 
since the magnetization vanishes identically at the crit¬ 
ical point with no background value (see Appendix [b|) . 
Without such terms, these approximants are also known 
as d-log Fade approximants. The order O now equals 
M + N. 

The boundary magnetization vanishes with the scaling 
dimension of the boundary operator (with iz = 1 ) 

= la <=“> 

Figure |^a) and Qb) show corresponding results calcu¬ 
lated up to order (/i/J)^° for the magnetization. Note, 
however, that the magnetization series only has even 
powers, so it is effectively a ten-term series in (/i/J)^. 

Notice that for the magnetization, the variation of the 
estimated exponents with order appears to adhere more 
closely to a dependence. Such terms are always 

present in addition to the terms, but should usually 
be weaker. However, we must use even more caution in 
extrapolating the magnetization exponents to asymptotic 
values, especially for a > 2 , as our data is still very far 
from the final values. 


IV. RENYI ENTROPIES IN D = 2 + 1 CORNERS 

Having seen the importance of unusual corrections to 
scaling for Renyi quantities in space-time dimension H-l, 
we turn to a case of intense recent interest: universal 
Renyi entropy corrections for D = 2 1 critical sys¬ 

tems. In two spatial dimensions, universal quantities ap¬ 
pearing in the Renyi entropies have a rich dependence 
on subregion geometry. For example, various universal 
terms (sub-leading to the area law) can arise from ge¬ 
ometries such as smooth circular bipartitions mm-, 
sharp corners [H4H57] ; bipartitioned infinite cylinders [10] 
or cylindrically-bifurcated tori [551l61j . 

Significant computational effort has been dedicated re¬ 
cently toward calculating the subleading contributions 
to the entanglement entropy arising from tt/2 corners in 
quantum models on the square lattice. Here we study 
the free scalar field theory, as well as the two-dimensional 
Ising chain in transverse field. A bipartition with a sin¬ 
gle corner (e.g. a quadrant of the infinite plane) in the 
entangling boundary contributes a subtractive, divergent 
logarithmic correction to the Renyi entropy, 

-S'a = - aa(7r/2) log ^0 H-, (32) 

where £ is the linear length scale of the bipartition A, and 
6 is the UV cutoff corresponding to the lattice length- 
scale. The constant Uq is not polluted by UV details, 
and constitutes a universal number that can be explored 
easily for tt/2 using lattice numerics. 

Previous numerical studies reveal that the value of 
aa(7r/2) obtained for the 0{N) model is somewhat close 
to N times the value obtained for a free scalar held the¬ 
ory by Casini and Huerta |54] . In this section, we use the 
numerical linked-cluster expansion (NLCE) to examine 
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FIG. 4. (a) Magnetization exponents estimated with the full length series. The A^-axis labels different approximates. The 
solid circles are the estimated values (with dashed lines as a guide to the eye). The solid lines are the theoretically predicted 
values, (b) Estimated exponents for different orders O of the series expansion. The squares are the asymptotic theoretical 
values. Dashed lines are guides to the eye, showing that a rough extrapolation of data is consistent with theory. 


the finite-size scaling behavior of Oo,(7r/2) for free bosons, 
and perform a parallel comparison of data obtained for 
the transverse field Ising model (TFIM) in two spatial 
dimensions, at its TV = 1 fixed-point. For the latter we 
use the density matrix renormalization group (DMRG) 
as our cluster solver [SH [SS]. 

In our NLCE procedure, described in detail in Ap¬ 
pendix the leading-order area-law piece of Eq. (32) 
cancels, leaving only the logarithmic divergence, 


Va = -aa \0g£ + ba. 


(33) 


In the next two sections, we interpret the order O of the 
NLCE (see Appendix as a length scale a procedure 
used successfully several times in the past This 


length scale is used to fit the two parameters of Eq. (33) 
to extract the universal constant ao,(7r/2). 


A. Free scalar field theory 


We begin by employing the NLCE method to study the 
free scalar field theory, described in the Appendices. For 
the free boson, the NLCE procedure can be used with a 
correlation-matrix technique as the cluster solver, as de¬ 
tailed in Appendix The correlation-matrix technique 
permits the calculation of Renyi entropies on very large 
finite-size clusters, allowing the NLCE sum to be carried 
to extremely high order. 

In Fig. [^a) , we extract the universal coefficient 
aa{'K/2) from NLCE data over a range of orders, from 


Omin to Omax, by fitting Eq. (331 as a function of NLCE 


order with two free parameters (oq, and ha)- For refer¬ 
ence, Fig. I^a) includes the exact values for the thermo¬ 
dynamic limit obtained by Casini and Huerta |54| . The 
von Neumann corner coefficient ai(7r/2) converges quite 


rapidly to the exact value, while the Renyi coefficients 
02 ( 71 / 2 ) and 03 ( 71 / 2 ) converge much more slowly, even 
for very large linear cluster sizes Omax = 50. 


As discussed in Section |ITj one may expect that for 
a > 1 a second extrapolation could become necessary be¬ 
cause of the unusual corrections to scaling arising from 
the conical singularities. We perform this second extrap¬ 
olation in Fig. [^b). The individual points are obtained 
by first fitting the NLCE data as a function of order to 
Eq. (331 from Omin up to a fixed maximum order Omax- 
This is the same procedure discussed above for produc¬ 
ing the curves labeled NLCE in Fig. [^a), but keeping a 
fixed and combining results for various Omax into a single 
curve. Next, we identify each Omax = L, a length scale 
which depends on the NLCE scheme (see Appendix [O. 
We then fit the finite-size results to the function A+B/L^ 
with A,B and 9 as fitting parameters. From Eq. (|^, 
taking a:£; = lforO = 2-l-l, one might guess the correc¬ 
tion should go as i.e. 9 = 1/a, in analogy to the 

D = 1-1-1 result. However, as we see in Fig.j^b), the opti¬ 
mized values of 9 which produce the best fit are somewhat 
different from 1 /a and yield extrapolations closer to the 
exact results [S3] , likely indicating different scaling expo¬ 
nents associated with the higher-dimensional system. 


From these extrapolations, we estimate the follow¬ 
ing values for aQ,(7r/2) in the thermodynamic limit: for 
a = 2 we obtain 02 ( 71 / 2 ) = 0.0065 which is 1% off of 
the exact value 0.0064 of Ref. |S3]; for a = 3 we obtain 
03 ( 71 / 2 ) = 0.0053 which is 3% from the exact value of 
0.0051. Clearly, our NLCE estimates will approach the 
exact results with increasing accuracy as Omax increases 
to infinity. 
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FIG. 5. (a) Numerical linked cluster results for the universal corner coefficient of the free scalar field theory, as a function of 
Renyi index a. Each curve was obtained by fitting NLCE data to the form Eq. (331, using only data between the minimum 
order and maximum order indicated, (b) Two Renyi coefficients, obtained by fitting NLCE results from order 3 to order L, as 
a function of maximum order L. Solid curves are fits to A + BjL^, with coefficients in the key obtained by fitting the largest 
10 L points. 


B. Transverse-field Ising model 

We turn now to NLCE calculations of the transverse- 
field Ising model (TFIM) in spatial dimension d = 2. The 
Hamiltonian on a square lattice is 

H = -/r J (34) 

and we work at the quantum critical point h/ J « 3.044. 
The corner contribution to the Renyi entropy has been 
studied by NLCE several times in the past [551EO] • The 
main technical difference from the free boson calculation 
of the preceding section is that each individual cluster 
must be solved using a method for strongly correlated 
Hamiltonians. This severely restricts the value of Omax, 
even if one uses the powerful DMRG method as a cluster 
solver. 

In Fig.|^ the universal coefficient aa{n/2) is extracted 
from NLCE data in a range of orders, from Omin to 
Omax, using the fitting form Eq. ( [3^ with two free pa¬ 
rameters (oq, and 6 q,). Up to the maximum order we 
reach, the results obtained this way for the TFIM are 
significantly different from the exact results for the free 
boson, particularly for a > 2. However, if one performs 
a second finite-size extrapolation for the Renyi coeffi¬ 
cients, using the form A + B/L^ with 0 values taken to be 
^TFiM = 1-41 X Ob ^ where 9b are the optimal exponents 
for the free boson and « 1.41 is the TFIM energy op¬ 
erator scaling dimension Eq. then our estimates for 
aa(7r/2) in the limit L —>■ oo are 0.0059 for a = 2; 0.0045 
for a = 3; and 0.0040 for a = 4. As shown in Fig. |^a), 
if we instead assume the exponents 6>tfim = 1.41/a, our 
estimates are 0.0061 for a = 2; 0.0048 for a = 3; and 
0.0045 for a = 4, higher than the extrapolation based on 


the boson exponents, especially for larger a. 

The data indicate that a large part of the observed dif¬ 
ference between free bosons and the TFIM in past stud¬ 
ies may be the unusual corrections to scaling, which are 
manifest as large finite-size effects in the NLCE proce¬ 
dure for a ^ 1. When such unusual corrections, arising 
from the presence of the conical singularity in the Renyi 
entropies, are taken into account, it becomes much more 
difficult to distinguish the value of aQ,(7r/2) obtained for 
the TFIM from that for free bosons. Since it will be very 
difficult to obtain larger values of O^ax using the DMRG 
method, it is hard to judge how much of the remaining 
discrepancy can be attributed to the very limited cluster 
sizes of the TFIM NLCE calculation. 


C. Discussion 

In this paper, we have performed a systematic study of 
the finite-size scaling of Renyi entropies and Renyi cor¬ 
relators at a quantum critical point. In particular, spe¬ 
cial attention was paid to the first subleading correction, 
which can play a crucial role in a proper extraction of 
universal terms. 

In space-time dimension I-I-1, a known scaling ansatz 
has been derived by Calabrese and Cardy and applied 
previously to numerical studies of Renyi entropies. We 
have applied this ansatz to Renyi correlators, calculated 
with conformal field theory, a mapping to free fermions, 
and series expansions for the transverse-field Ising model 
(TFIM) at its quantum critical point. For data obtained 
with finite-size (or finite-series length) calculations, we 
find that a second extrapolation in the size (or length 
of the series) using this ansatz is needed to show con- 
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FIG. 6. (a) Numerical linked cluster results for the universal corner coefficient at the quantum critical point of the TFIM, 
as a function of Renyi index a. The lower curves were obtained by fitting NLCE data to the form Eq. (331, using only data 
between the minimum order and maximum order indicated. The two curves labeled “Extrapolation” were obtained using the 
two-step fitting procedure described in the text, in one case hxing the power-law exponents to be 0 = 1.41/a and in the other 
case 6 = 1.41 • 6b where 6b is the optimal exponent for fitting the finite-size free boson results, (b) Finite-size Renyi corner 
coefficients estimated fitting NLCE data up to some maximum order L. Solid curves are fits to A + B/L^ with 6 — 1.41 ■6b- 
The intercepts at l/L —> 0 define the solid curve labeled “Extrapolation” in subfigure (a). 


sistency with the theoretically known values. Failing to 
take the correction to scaling into account can lead to 
qualitatively incorrect answers. Thus, it appears that the 
conical singularities are far more important for Renyi cor¬ 
relators in D = 1-1-1 than irrelevant operators are, e.g. for 
computing bulk properties of the 3-dimensional classical 
Ising model. 

We also studied 2-1-1 space-time dimensions, where the 
influence on scaling of the conical singularity is generally 
not known exactly. We focus on the universal term Ua 
which arises in Renyi entropies due to a 7r/2 corner, using 
a recently developed numerical linked cluster expansion 
(NLCE). In a scalar field-theory corresponding to free 
bosons, where the exact result is known, very large sys¬ 
tems can be studied. As in 1-1-1, we see that strong cor¬ 
rections to scaling are present for a > 1 Renyi entropies, 
which retain a significant discrepancy from exact results 
even for relatively large cluster sizes. However, if we do 
an additional extrapolation in the maximum linear di¬ 
mension as L~^, we find that fitting 0 as a function of a 
produces results that agree well with the exact value. 

Finally, following a similar procedure for the corner 
Renyi entropies of the 2-1-1 TFIM at its quantum critical 
point, we extrapolate existing NLCE data. Using the 
optimal 0 found for free bosons, but including a non¬ 
trivial scaling dimension for the energy operator, we find 
that the values for aa{Tr/2) are much closer to the free 
boson values than concluded from previous analyses. 

Improving numerical accuracy with such extrapola¬ 
tions is crucially important for the interpretation of many 
lattice calculations of Renyi entropies. For example, re¬ 
cent calculations on several models tuned to the 0{N) 
Wilson-Fisher fixed-point have revealed that aQ,(7r/2) is a 


universal function of a that scales to within numerical ac¬ 
curacy as aa(7r/2) ^ Nca{TT/2), where N is the number of 
field components, and Cq is some universal function that 
appears to be the same for all 0{N) models studied to 
date. Intriguingly, the central charge Ct also grows lin¬ 
early with N up to small corrections [131 El] ■ Our results 
indicate that, when unusual corrections to scaling are 
taken into account, 00 ( 71 / 2 ) for = 1 can be very close 
numerically to the value for free bosons, however for other 
0{N) values, the conclusion that aa{Tr/2) ^ Nca{n/2) 
remains well justified. Thus, we confirm that the univer¬ 
sal subleading corner coefficient is sensitive to number of 
degrees of freedom in the low-energy field theory describ¬ 
ing a strongly-interacting quantum system at its critical 
point. 
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Appendix A: Reduced density matrices and 
entropies for free theories 


We recall in this appendix some general results re¬ 
garding reduced density matrices for free bosonic and 
fermionic systems. The techniques are well-known; we 
refer to Ref. [M] for a review. We first present the 
method to compute entanglement for free bosons, rele¬ 
vant to Sec. IV A before focusing on Renyi correlations 
(relevant to Sec. IIIB). 


1. Entropies for the free lattice scalar field 


Since (O) = Tr (paO) for any observable O in subsystem 
A, the reduced density matrix itself is the exponential of 
a quadratic boson form 



<^gf>lbq 


(AS) 


where the operators bq, are linear combinations of the 
original fields (j)i,TTi. They also satisfy the canonical com¬ 
mutation relations [bq,b^q,] = Sqqi. Here Z is a. normaliza¬ 
tion constant that ensures Tr(pyi) = 1. Such an ansatz 
has by definition the Wick factorization property build 
in. The next step is then to adjust the energies e/ and 
the operators bj so as to reproduce the ground-state cor¬ 
relations, which can always be done. After some algebra, 
we find that the single particle energies ei satisfy 


It has been long known how to evaluate the entangle¬ 
ment entropy for a system of coupled oscillators, see e.g. 
Refs. [65l[66]- Here we follow the method put forward 
by Peschel [^, that was used by Casini and Huerta [S3] 
to compute exactly the corner contribution mentioned in 
the text. Numerically, the technique may also be used to 
compute the entropy on very large clusters. 

Consider a finite two-dimensional square lattice, such 
that a free scalar field (j)i and its conjugate momentum 
TTi exist at each lattice point. A non-interacting system 
is described by following lattice Hamiltonian, 


^ — n ^ ( ^X,y T i4‘x+l,y 4’x,y) + {(j^x^y+l 'Px,y) 


x,y—l,l 


+ rn^4>l,y 


(Al) 


where and Ny are the linear dimensions of the lattice. 
The Hamiltonian can be rewritten in terms of bosonic 
operators, but we refrain from doing so for now. We are 
mainly interested in the massless case m = 0 , but keep 
general m in all the formulas. 

The total number of sites is N^Ny. The Hamiltonian 
may be rewritten as 

N^Ny N 

H = - (jiiMijcf), (A2) 

2=1 


(l/ 2 )coth(eg/ 2 ) = y, 




(A 6 ) 


where the Vq are the positive eigenvalues of the matrix 
Ca = VXaPa- Here Xa (resp. Pa) is the matrix X 
(resp. P) restricted to subsystem A. We note it Ca, to 
emphasize the fact that it depends on the choice of sub¬ 


system. The Renyi entropy then follows from Eq. (A5) 


= 


1 


Of — 1 


Tr log [{Ca + 1/2)“ - {Ca - 1/2)“]. (A7) 


This is of course a huge simplification, because what 
is only needed is the diagonalization of the matrix Ca, 
whose size is given by the number of sites in A. This last 
step can be done using standard linear algebra routines. 

For the geometry we choose the full system to be an 
open-boundary rectangle of size x Ny. The correla¬ 
tors come from diagonalizing the discrete laplacian on 
the rectangle, and using Eqs. (A3 A4). We obtain 


{(t>x,y(l)x-,y') = (jv +1){N -H) ^ sm{q,,x)sm{qq,x') 

(A8) 

^{qx,qy) 

T^x,yT^x',y') = sin{q,:x)sm{q,,x') 

' ^ Ay ' qq,,qy 

X sin {qyy) sin (%?/') x ui{qx,qy) (A9) 


X sin {qyy) sin(gy?/') 
2 


The matrix M is, up to a diagonal term, nothing but 
the discrete laplacian. The ground state correlations are 
given 


Xij — {4>i4>j) — 2 (A^ ^ )b’ 

(A3) 

Pij = = ^(A^ ^ 

(A4) 


and may be obtained explicitly in many geometries, nu¬ 
merically in others. 

Let us now cut our system in two parts A and B, and 
focus only on the degrees of freedom inside subsystem A. 


where the quasi-momenta kx, ky are quantized as 

J2, 'JY 

= ’ nx = l,2,...,Nx, (AlO) 

TL TT 

~ ]\l _|_ 2 ^ ’ ny = 1,2,Ny. (All) 

and 

uj{qx,qy) = ^j4sin^{qx/2) + 4sm‘^{qy/2) +m'^. (A12) 

Combined with the NLCE procedure described in Ap¬ 
pendix this produces the data shown in Sec. |IV A| 
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We note that the choice of boundary conditions is not 
completely innocent. For example, if we were to choose 
a periodic-boundary torus instead of an open-boundary 
rectangle, the correlators would be given by similar but 
translation invariant formulas 


{(t>0,o4>x,y) 9 AT /V 


COS (kxx) COS (kyu) 


(A13) 


2N^Ny ^ Lo{koo,ky) 

{'^0,0'^X,y) = yy- ■E cos (kxx) cos {kyu) oj{kx, ky), 

(AM) 


kx,ky 


with 


The coefficients akj , Pkj and the single particle energies 
Cfc may be e.g. obtained by diagonalizing the following 
2L X 2L Bogoliubov matrix 


M = 


U V \ 
V* 1-U* ) 


(A21) 


where 17* (resp. V*) is the transpose of U (resp. con¬ 
jugate transpose of V). Using this method, computing 
for example the ground-state correlations is a straight¬ 
forward task. 

We now focus on a subsystem, which we choose to be 
the first € sites for simplicity. For later convenience we 
introduce the following 2£ x 2^ correlation matrix 


uj{kx, ky) = \J4:SiT?(kx/2) + As\T?{ky/2) w?. (A15) 


K = 


C D \ 
D* 1-C* J ' 


(A22) 


and the momenta quantized as 


2nxTr 

Nx 


.,Nx-l, (A16) 

2nyTT 

Ny 

^ 5 ■ 

.,Ny-l. (A17) 

Unfortunately, the correlator (AfS) 

is divergent in the 


massless case m = 0, due to the presence of a zero mode 
{kx,ky) = (0,0). Thus, Sa also diverges, as shown pre¬ 
viously in Ref. [2], significantly complicating the calcu¬ 
lation of Renyi entropies on finite-size tori. To remove 
this divergence, a possible way is keep a very small mass 
term; this results in an additional term proportional to 
log (mL) in the entropy, where L is the length of the 
boundary. This m may then be systematically decreased 
and the results extrapolated to to —>■ 0 . 

This complication is not a fundamental one; however, 
it motivates us to only consider clusters with open bound¬ 
ary conditions in the procedure of Appendix which do 
not have zero modes. 


2. Renyi correlations 


C and D encode the two types of correlators: 

C = {{clcj))i<ij<e, (A23) 

D = ((clcj.))i<,,,<,. (A24) 

Now, any correlation in subsystem A may be recovered 
from the reduced density matrix, 

{0)=Ty{paO). (A25) 


Because of Wick’s theorem, the reduced density matrix 
itself can be written as 


PA 


Tr ’ 


(A26) 


where Ha is a quadratic fermion form for a system of 
size i, similar to (A18). We name Ha the entanglement 
Hamiltonian, as it is different from the original Hamilto¬ 
nian H. The crucial point is that the Bogoliubov matrix 
corresponding to Ha and the subsystem correlation ma¬ 
trix are diagonalized by the same transformation [48j . 
More precisely, it is possible to find a unitary matrix U 
which diagonalizes K 


We consider the following Hamiltonian 

L 

77 = ^ (Uijclcj + -I- h.cj . (A18) 

*j=i 

Such a quadratic fermion form can always be diagonal¬ 
ized by a Bogoliubov transformation. Namely, there is a 
set of fermion operators 


_ rr* f ^^agg (Ag) 0 

“ V 0 diag, (1 - A,) 


U, 


and such that 




9=1 


with 


(A27) 


(A28) 


3 

which obey the canonical commutation relations 
{4)'7fc'} = 5kkA and which diagonalize 77: 


= U(cl,...,cj,ci,...,cf)*. 

(A29) 

The single particle eigenenergies Sq of the entanglement 
Hamiltonian may then be obtained from the \q by re¬ 
quiring consistency with (A25). We obtain 


77 = ^ ekdldk (A20) 

k 


\ 


1 

f -I- e^i 


(A30) 
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This relation, together with (A29), uniquely determines 


the entanglement Hamiltonian from the eigenvalues and 
the eigenvectors of the correlation matrix K. 

With this correspondence at hand, computing Renyi 
correlations in subsystem A becomes straightforward. In¬ 
deed the new (powered up) reduced density is 


(a) 

Pa = 


{paY 


c—oHa 


Tr [{pa)Y Tr [e 


(A31) 


Therefore, the new entanglement Hamiltonian is simply 
uHa ■ We now wish to obtain the new correlation matrix 




Ca \ 


with 


Ca — 

= ((c|c|)^)l<i , 

corresponding to the Renyi correlators 


(O) 


Tr 




(A32) 


(A33) 

(A34) 


(A35) 


To do so it is sufficient to reverse the logic leading to the 
determination of Ha- From (A27l and (A30), we have 


K = U* 


f diag, (t^) 

\ 0 diag, ( 


0 


l+e- 


U (A36) 


so the new correlation matrix may be obtained by 
simply multiplying all single particle energies Sq by a. 
We obtain 


K^ = U* 



) 



(A37) 

In terms of the initial correlation matrix AT, the previous 
relation reads 


K^=[l + (K-^ - I)“] 


-1 


(A38) 


The result is true for any eigenstate and also at finite 
temperature, the only requirement being a free fermion 
Hamiltonian. In fact, (A38) even holds for quadratic 


boson Hamiltonians, provided the correlation matrix be 
redefined as 


K = 


C* D \ 
D* 1 + C* ) ' 


(A39) 


instead of ( |A22 ). This can be shown by repeating the 
steps described above with bosonic operators instead of 
fermions. 


Appendix B: Series analysis 

Series are analyzed by the method of Differential Ap- 
proximants These are generalizations of the bet¬ 

ter known d-log Fade method and allow for a power-law 
singularity. The function of interest f{x) is represented 
by a solution to a first order ordinary differential equa¬ 
tion: 


= Pi{x) (HI) 

Here, P^{x), P^{x) and P^(x) are polynomials of order 
M, N and J respectively that are determined by demand¬ 
ing that the series expansion for the function f{x) match 
to some order. It is easy to show that these differential 
equations have power-law singularities at critical points 
given by the roots of the polynomial P^ {x) with expo¬ 
nents related to the residue of Pf^ (x) and the value of 
PY {x) at that critical point. The less important P^(x) 
term allows for a background function. Since the critical 
point of Eq. ([^ is known to be at J//i = 1, we further 
bias the analysis by demanding that 

P,^{x = l)=0, (B2) 

thus ensuring that critical point is at ccc = 1. 

The calculation proceeds by first choosing orders for 
the polynomials M, N and J. We fix J = 0 that just 
allows for a constant background term, and keep M to 
at least 2. The sum of the orders {M + N -\- J) can not 
exceed the number of terms (more correctly number of 
terms minus one after the biasing) that one knows in the 
series expansion. In general, different choice of the order 
of the three polynomials M, N and J will give different 
values for critical exponents. 


Appendix C: The Numerical Linked Cluster 
Expansion 

The procedure of the Appendix [fallows one to calcu¬ 
late the entanglement entropy for free bosons on clusters 
of arbitrary shape and size. Thus, it may be used as a 
“cluster solver” for a NLCE procedure, similar to previ¬ 
ous studies where Lanczos or DMRG were used as cluster 
solvers for interacting models p5H5T| . 

We refer the reader to the relevant literature [BTUTCI] . 
The NLCE method is based on the fact that an exten¬ 
sive property P of a lattice model can be expressed as 
a sum over contributions from all distinct clusters which 
are embeddable in the lattice. Then, this property per 
site is, 


P/iV = ^W(c), (Cl) 

c 


where W (c) is the “weight”, or the unique contribution 
to P that comes from a cluster c m- It is defined recur- 


















13 


sively, 


in Eq. (32)). As described in Ref. [53] the property P 


W{c) = P{c) 


(C2) 


sGc 


where s is any subcluster of c. As dictated by the 
inclusion-exclusion principle, only connected clusters give 
non-zero weights to the sum. The NLCE procedure 
builds up the value of P/N starting from the smallest 
cluster c = 1 and ending with some maximal cluster size. 
For this paper, we employ only rectangular clusters [53] . 
which makes the calculation of cluster embeddings s € c 
trivial; then, the computational bottleneck is shifted to 
the cluster solver. 

For interacting models, the computational cost of cal¬ 
culating the given property P on the cluster is extremely 
expensive. However, as described in the last section, a 
very efficient numerical method can be devised for calcu¬ 
lating the Renyi entropies for a free bosons on the square 
lattice using Equations ( A8|A9 ). This allows us to per¬ 
form the NLCE to extremely high orders. 

As with any cluster solver on the square lattice, the 
NLCE is able to isolate the entanglement due to a 7r/2 
corner by subtracting off the contributions to Sa from 
the linear portions of the boundary (the area-law term 


then becomes the isolated contribution for a single cor¬ 
ner, given by Eq. (33). This property is made extensive 


by adding the values for all possible translations of the 
corner around the cluster c. 

Finally, for studies of critical systems such as the free 
boson, the definition of a length scale is important. In 
this paper, we are interested in benchmarking the success 
of the NLCE in converging the exact results for aQ,(7r/2) 
in the limit where this length scale approaches infinity. A 
second advantage of using rectangular clusters is the easy 
definition of this length scale in relation to the cluster or¬ 
der O. In past studies, several definitions have been used, 
including the arithmetic and geometric means of the rect¬ 
angle linear dimensions IMMlIZa- In the main text, we 
define the cluster length scale ^ as the maximum of the 
two linear dimensions. This definition, which also defines 
the expansion order O, is simple to use and is motivated 
by the intuition that the role of the smaller clusters in 
the NLCE sum is to cancel boundary effects. Thus long 
and narrow clusters of dimension N^^Ny should only be 
included at the same order as N x N square clusters with 
N = max{Nx, Ny); otherwise these narrow clusters are 
not the boundary of any more isotropic cluster appearing 
to this order in the sum. 
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